********************************************************************************************************************************
***   Replication file for:                                                                                                  ***
***   Berbee, P., Braun, S. T. and Franke, R. (2024). Reversing Fortunes of German Regions, 1926-2019. JOEG				     ***
***   							                                                                                             ***
***   SCRIPT: 	4_lopsided_economic_structure.do																			 ***	
***   PURPOSE: 	Generates the figures and tables in Section 4 ("Lopsided economic structure and limited adjustment			 *** 
***				capacity") and the associated Appendix A-3.																	 ***
***	  Tables:	2, 3, A-5, A-6, A-7, A-8										 										     ***
***	  Figures:	6, A-3, A-4, A-5																							 ***
********************************************************************************************************************************


* Preamble (unnecessary when executing run.do)
run "$reversing/scripts/programs/_config.do"

************
* Code begins
************

use "$reversing/processed/workingdataset.dta", clear

********************************************************************************
*** Table 2: Early industrialization and industry structure in 1907 (2SLS estimates)

local covariates log_land_access1 town_1700_perarea

** Column 1: Average firm size
* Conley
eststo twosls_firmsize_co: acreg firmsize_industry_1907 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
* Cluster
eststo twosls_firmsize_cl: ivreg2 firmsize_industry_1907 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)


** Columns (2)-(4): Employment in sectors with at least 50% in large firms; HHI index 
foreach outcome in f501 f1000 hhi {
* Conley
eststo twosls_`outcome'_co: acreg empshare_`outcome'_ind_1907 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
* Cluster
eststo twosls_`outcome'_cl: ivreg2 empshare_`outcome'_ind_1907 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)
}

** Column (5): Principal component
corr firmsize_industry_1907 empshare_f501_ind_1907 empshare_f1000_ind_1907 empshare_hhi_ind_1907
pca firmsize_industry_1907 empshare_f501_ind_1907 empshare_f1000_ind_1907 empshare_hhi_ind_1907
predict pcomp*

** 2SLS
* Conley
eststo twosls_pcomp1_co: acreg pcomp1 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
* Cluster
eststo twosls_pcomp1_cl: ivreg2 pcomp1 (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)

** Mediation analysis
ivmediate change_perc_5719 log_land_access1 town_1700_perarea if year == 2019, mediator(pcomp1) treatment(empshare_ind_1882_std) instrument(log_coal_access1) vce(cluster rb_id) full

** Export to Latex
esttab twosls_firmsize_co twosls_f501_co twosls_f1000_co twosls_hhi_co twosls_pcomp1_co using"$reversing/results/tables/tab2_part1.tex", keep (empshare_ind_1882_std) coeflabels(empshare_ind_1882_std "Employment share industry 1882") cells(b(star fmt(3)) se(fmt(3) par("(" ")"))) nomtitles nonum noobs collabel(none)  star(* 0.10 ** 0.05 *** 0.01) plain fragment booktabs replace

esttab twosls_firmsize_cl twosls_f501_cl twosls_f1000_cl twosls_hhi_cl twosls_pcomp1_cl using"$reversing/results/tables/tab2_part1.tex", keep(empshare_ind_1882_std) coeflabels(empshare_ind_1882_std " ") cells(se(fmt(3) par("[" "]")))  nomtitles nonum noobs collabel(none) plain fragment booktabs append	

*** Means and standard deviations
estpost tabstat firmsize_industry_1907 empshare_f501_ind_1907 empshare_f1000_ind_1907 empshare_hhi_ind_1907 pcomp1 if year == 2019, c(stat) stat(mean sd)
*ereturn list
matrix desc = e(mean) \ e(sd)
matrix list desc
esttab matrix(desc, fmt(3)) using"$reversing/results/tables/tab2_part2.tex", coeflabels(mean "Mean" sd "Standard deviation") nomtitles collabel(none) plain  booktabs fragment replace
matrix drop desc



********************************************************************************
*** Table 3: Early industrialization and adjustment pressure (2SLS estimates)


local covariates log_land_access1 town_1700_perarea

** 2SLS

* Clustered standard errors
eststo ind_cl: ivreg2 empshare_ind_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo montan_cl: ivreg2 empshare_montan_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo modern_cl: ivreg2 empshare_modernindustries_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo other_cl: ivreg2 empshare_otherindustry_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo ind_2019_cl: ivreg2 empshare_prod_industry `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)

* Conley standard errors
qui eststo ind_co: acreg empshare_ind_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
qui eststo montan_co: acreg empshare_montan_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
qui eststo modern_co: acreg empshare_modernindustries_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
qui eststo other_co: acreg empshare_otherindustry_1950 `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
qui eststo ind_2019_co: acreg empshare_prod_industry `covariates' (empshare_ind_1882_std=log_coal_access1) if year== 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)

** Export to latex

esttab ind_co montan_co modern_co other_co ind_2019_co using"$reversing/results/tables/tab3_part1.tex", keep (empshare_ind_1882_std) coeflabels(empshare_ind_1882_std "Employment share industry 1882") cells(b(star fmt(3)) se(fmt(3) par("(" ")"))) nomtitles nonum noobs collabel(none)  star(* 0.10 ** 0.05 *** 0.01) plain fragment booktabs replace

esttab ind_cl montan_cl modern_cl other_cl ind_2019_cl using"$reversing/results/tables/tab3_part1.tex", keep(empshare_ind_1882_std) coeflabels(empshare_ind_1882_std " ") cells(se(fmt(3) par("[" "]")))  nomtitles nonum noobs collabel(none) plain fragment booktabs append	


*** Means and standard deviations
estpost tabstat empshare_ind_1950 empshare_montan_1950 empshare_modernindustries_1950 empshare_otherindustry_1950 empshare_prod_industry if year == 2019, c(stat) stat(mean sd)
*ereturn list
matrix desc = e(mean) \ e(sd)
matrix list desc
esttab matrix(desc, fmt(3)) using"$reversing/results/tables/tab3_part2.tex", coeflabels(mean "Mean" sd "Standard deviation") nomtitles collabel(none) plain  booktabs fragment replace
matrix drop desc


********************************************************************************
* Table A-5: 2SLS mediation analysis

* Panel A: Firm size 1907
ivmediate change_perc_5719 log_land_access1 town_1700_perarea if year == 2019, mediator(firmsize_industry_1907) treatment(empshare_ind_1882_std) instrument(log_coal_access1) vce(cluster rb_id) full

* Panel B: Large firms > 501 workers 1907
ivmediate change_perc_5719 log_land_access1 town_1700_perarea if year == 2019, mediator(empshare_f501_ind_1907) treatment(empshare_ind_1882_std) instrument(log_coal_access1) vce(cluster rb_id) full

* Panel C: Large firms > 1000 workers 1907
ivmediate change_perc_5719 log_land_access1 town_1700_perarea if year == 2019, mediator(empshare_f1000_ind_1907) treatment(empshare_ind_1882_std) instrument(log_coal_access1) vce(cluster rb_id) full

* Panel D: HHI index of industry concentration 1907
ivmediate change_perc_5719 log_land_access1 town_1700_perarea if year == 2019, mediator(empshare_hhi_ind_1907) treatment(empshare_ind_1882_std) instrument(log_coal_access1) vce(cluster rb_id) full




********************************************************************************
* Table A-6: Effects of early industrialization on the population share of enrolled students, 2019

eststo clear

local covariates log_land_access1 town_1700_perarea
local outcomelist stud_total_pop_2019 stud_humanities_pop_2019 stud_law_econ_pop_2019 stud_math_natsc_pop_2019 stud_medicine_pop_2019 stud_engineer_pop_2019 stud_other_pop_2019

foreach outcome in `outcomelist'{
eststo `outcome'_cl: ivreg2 `outcome' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo `outcome'_co: acreg `outcome' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
}


** Export results to Latex
esttab stud_total_pop_2019_co stud_humanities_pop_2019_co stud_law_econ_pop_2019_co stud_math_natsc_pop_2019_co stud_medicine_pop_2019_co stud_engineer_pop_2019_co stud_other_pop_2019_co using"$reversing/results/tables/tabA6_part1.tex", keep (empshare_ind_1882_std) star(* 0.10 ** 0.05 *** 0.01) coeflabels(empshare_ind_1882_std "Employment share industry 1882") se(%9.3fc) b(%9.3fc) parentheses replace nonotes nodepvars fragment plain noobs nomtitles nonum collabels(none) booktabs
esttab stud_total_pop_2019_cl stud_humanities_pop_2019_cl stud_law_econ_pop_2019_cl stud_math_natsc_pop_2019_cl stud_medicine_pop_2019_cl stud_engineer_pop_2019_cl stud_other_pop_2019_cl using"$reversing/results/tables/tabA6_part1.tex", keep(empshare_ind_1882_std) coeflabels(empshare_ind_1882_std " ") cells(se(fmt(3) par("[" "]"))) prefoot("\midrule") nomtitles nonum noobs collabel(none) plain fragment booktabs append

* Means and standard deviations
estpost tabstat stud_total_pop_2019 stud_humanities_pop_2019 stud_law_econ_pop_2019 stud_math_natsc_pop_2019 stud_medicine_pop_2019 stud_engineer_pop_2019 stud_other_pop_2019  if year == 2019, c(stat) stat(mean sd)
matrix desc = e(mean) \ e(sd)
matrix list desc
esttab matrix(desc, fmt(3)) using"$reversing/results/tables/tabA6_part2.tex", coeflabels(mean "Mean" sd "Standard deviation") nomtitles collabel(none) plain  booktabs fragment replace
matrix drop desc




********************************************************************************
* Table A-7: Effects of early industrialization on educational attainment
* University, at least higher secondary or vocational, no school leaving qualification

eststo clear
rename popshare_noschooldegree* share_nosdegree*

local covariates log_land_access1 town_1700_perarea
local outcomelist share_high_ed share_med_high_ed share_nosdegree

foreach year in 1970 2011 {
foreach outcome in `outcomelist'{
eststo `outcome'_`year'_cl: ivreg2 `outcome'_`year' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
eststo `outcome'_`year'_co: acreg `outcome'_`year' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
}
}

** Export results to Latex
esttab share_high_ed_1970_co share_high_ed_2011_co share_med_high_ed_1970_co share_med_high_ed_2011_co share_nosdegree_1970_co share_nosdegree_2011_co using"$reversing/results/tables/tabA7_part1.tex", keep (empshare_ind_1882_std) star(* 0.10 ** 0.05 *** 0.01) coeflabels(empshare_ind_1882_std "Employment share industry 1882") se(%9.3fc) b(%9.3fc) parentheses replace nonotes nodepvars fragment plain noobs nomtitles nonum collabels(none) booktabs
esttab share_high_ed_1970_cl share_med_high_ed_1970_cl share_high_ed_2011_cl share_med_high_ed_2011_cl share_nosdegree_1970_cl share_nosdegree_2011_cl using"$reversing/results/tables/tabA7_part1.tex", keep(empshare_ind_1882_std) coeflabels(empshare_ind_1882_std " ") cells(se(fmt(3) par("[" "]"))) prefoot("\midrule") nomtitles nonum noobs collabel(none) plain fragment booktabs append

* Means and standard deviations
estpost tabstat share_high_ed_1970 share_high_ed_2011 share_med_high_ed_1970  share_med_high_ed_2011 share_nosdegree_1970 share_nosdegree_2011  if year == 2019, c(stat) stat(mean sd)
matrix desc = e(mean) \ e(sd)
matrix list desc
esttab matrix(desc, fmt(3)) using"$reversing/results/tables/tabA7_part2.tex", coeflabels(mean "Mean" sd "Standard deviation") nomtitles collabel(none) plain  booktabs fragment replace
matrix drop desc


********************************************************************************
*** Table A-8: Early industrialization, self-employment, and political outcomes

** Regressions
local covariates log_land_access1 town_1700_perarea
local outcomelist empshare_appr_ind_1970 hhi_industryemp_1950 empshare_self_1950 ob_years_dominant_party ob_years_SPD voteshare_SPD_1957
foreach outcome in `outcomelist'{
qui eststo `outcome'_co: acreg `outcome' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
qui eststo `outcome'_cl: ivreg2 `outcome' `covariates' (empshare_ind_1882_std=log_coal_access1) if year == 2019, cluster(rb_id)
}

** Export results to Latex
esttab empshare_appr_ind_1970_co hhi_industryemp_1950_co empshare_self_1950_co ob_years_dominant_party_co ob_years_SPD_co voteshare_SPD_1957_co using"$reversing/results/tables/tabA8_part1.tex", keep (empshare_ind_1882_std) star(* 0.10 ** 0.05 *** 0.01) coeflabels(empshare_ind_1882_std "Employment share industry 1882") se(%9.3fc) b(%9.3fc) parentheses replace nonotes nodepvars fragment plain noobs nomtitles nonum collabels(none) booktabs
esttab empshare_appr_ind_1970_cl hhi_industryemp_1950_cl empshare_self_1950_cl ob_years_dominant_party_cl ob_years_SPD_cl voteshare_SPD_1957_cl using"$reversing/results/tables/tabA8_part1.tex", keep(empshare_ind_1882_std) coeflabels(empshare_ind_1882_std " ") cells(se(fmt(3) par("[" "]"))) prefoot("\midrule") nomtitles nonum noobs collabel(none) plain fragment booktabs append


* Means and standard deviations
estpost tabstat `outcomelist' if year == 2019, c(stat) stat(mean sd)
matrix desc = e(mean) \ e(sd)
matrix list desc
esttab matrix(desc, fmt(3)) using"$reversing/results/tables/tabA8_part2.tex", coeflabels(mean "Mean" sd "Standard deviation") nomtitles collabel(none) plain  booktabs fragment replace
matrix drop desc



********************************************************************************
*** Figure 6: Impact of early industrialization on the rank in the patent per capita distribution, 1882-2016

*** Generate percentile ranks to ease interpretation
foreach period in 1878_1882 1903_1907 1922_1926 1931_1935 1951_1955 1962_1966 1970_1974 1983_1987 1992_1996 2002_2006 2012_2016 {
bysort year: egen patent_ranking_`period'=rank(pc_patents_`period')
gen patent_perc_`period' = round(patent_ranking_`period' / 1.63) // transform to percentile ranks
drop patent_ranking_`period' patent_ranking_`period'
}


* For plotting: transform periods in variable names to years
rename *1878_1882 *1882
rename *1903_1907 *1907
rename *1922_1926 *1926
rename *1931_1935 *1935
rename *1951_1955 *1955
rename *1962_1966 *1966
rename *1970_1974 *1974
rename *1983_1987 *1987
rename *1992_1996 *1996
rename *2002_2006 *2006
rename *2012_2016 *2016


local covariates log_land_access1 town_1700_perarea 
local years 1882 1907 1926 1935 1955 1966 1974 1987 1996 2006 2016

matrix C90 = J(3,11,.)
matrix C95 = J(3,11,.)
local j = 1
foreach year in `years' {
		acreg patent_perc_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100)
		lincom empshare_ind_1882_std, level(90)
		matrix C90[1,`j']= r(estimate) \ r(lb) \ r(ub)
		lincom empshare_ind_1882_std, level(95)
		matrix C95[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}


local colnames
foreach year in `years' {
    local colnames `colnames' `year'
}

matrix colnames C90 = `colnames'
matrix list C90

matrix colnames C95 = `colnames'
matrix list C95

coefplot (matrix(C95), ci((C95[2] C95[3])) mcolor(black) ciopts(lpattern(solid) lcolor(black) msize(medlarge) recast(rcap))) (matrix(C90), ci((C90[2] C90[3])) mcolor(black) ciopts(lpattern(solid) lcolor(gray) lwidth(medium) mcolor(gray) msize(medlarge) recast(rcap))), xtitle("Year") yline(0, lcolor(red)) vertical graphregion(color(white)) at(_coef) xlabel(1880(20)2020) ciopts(recast(rcap)) name(patentshare, replace) ytitle("Percentile ranks") legend(off)

graph export "$reversing/results/figures/fig6.eps", replace




********************************************************************************
*** Figure A3: The effect of early industrialization on the employment shares in industry and services, 1939-2019

* Note: In contrast to the other regression figures, we only report 90% confidence intervals as otherwise the figure may
* be difficult to read as it reports two sets of regression estimates.

* Note: Data sources are the GESIS files 1950-1987 as well as the GDP data. The latter contains sectoral employment data beginning in 1992. Here we choose 1994 instead, as the data on 1992 come from the 1980-1992 GDP revision.

local covariates = "log_land_access1 town_1700_perarea" 

matrix C = J(6,9,.)
local j = 1

local years "1939"
foreach year in `years' {
		ivreg2 empshare_ind_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)
		lincom empshare_ind_1882_std, level(90)
		matrix C[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}
local years "1950 1961 1970 1987"
foreach year in `years' {
		ivreg2 empshare_indust_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)
		lincom empshare_ind_1882_std, level(90)
		matrix C[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}
local years "1994 2000 2010 2019"
foreach year in `years' {
		ivreg2 empshare_prod_industry (empshare_ind_1882_std = log_coal_access1) `covariates' if year == `year', cluster(rb_id)
		lincom empshare_ind_1882_std, level(90)
		matrix C[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}

local colnames
local years "1939 1950 1961 1970 1987 1994 2000 2010 2019"
foreach year in `years' {
    local colnames `colnames' `year'
}
matrix colnames C = `colnames'
matrix list C


local years "1939 1950 1961 1970 1987"
local j = 1
foreach year in `years' {
		ivreg2 empshare_serv_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, cluster(rb_id)
		lincom empshare_ind_1882_std, level(90)
		matrix C[4,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}
local years "1994 2000 2010 2019"
foreach year in `years' {
		ivreg2 empshare_serv_total (empshare_ind_1882_std = log_coal_access1) `covariates' if year == `year', cluster(rb_id)
		lincom empshare_ind_1882_std, level(90)
		matrix C[4,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}

local colnames
local years "1939 1950 1961 1970 1987 1994 2000 2010 2019"
foreach year in `years' {
    local colnames `colnames' `year'
}
matrix colnames C = `colnames'
matrix list C

coefplot (matrix(C), ci((C[2] C[3])) ciopts(lcolor(gs4) recast(rcap)) msymbol(O) mcolor(gs4) label("Industry") ) (matrix(C[4]), ci((C[5] C[6])) ciopts(lcolor(gs12) recast(rcap)) label("Services") msymbol(T) mcolor(gs12)), xtitle("Year") ytitle("Employment share") yline(0, lcolor(red)) vertical graphregion(color(white)) at(_coef) xlabel(1940(10)2020) ciopts(recast(rcap))

graph export "$reversing/results/figures/figA3.eps", replace


********************************************************************************
*** Figure A-4: Mean industrial employment share by 1882 industrial employment, 1882-2019

preserve

summarize empshare_ind_1882 if year == 1926, d
gen empshare_ind_1882_m = inrange(empshare_ind_1882, 0, r(p50)) 
replace empshare_ind_1882_m = 2 if inrange(empshare_ind_1882, r(p50), 1) 

collapse (mean) empshare_ind_1882 empshare_ind_1907 empshare_ind_1939 empshare_ind_1950 empshare_indust_1950 empshare_indust_1961 empshare_indust_1970 empshare_indust_1987 empshare_prod_industry, by(empshare_ind_1882_m year)
drop if empshare_prod_industry==. & year<1980
drop empshare_ind_1950
rename empshare_indust_* empshare_ind_* 

reshape long empshare_ind_, j(jahr) i(empshare_ind_1882_m year)
replace year=jahr if year==1980
replace empshare_prod_industry=empshare_ind_ if year<1992
drop if year>=1992 & jahr>1882 /* keep just one obs per year */
drop empshare_ind_ jahr
drop if empshare_ind_1882_m==.

reshape wide empshare_prod_industry, i(year) j(empshare_ind_1882_m)

twoway (line empshare_prod_industry1 year, lcolor(gs4)) ///
 (line empshare_prod_industry2 year,  lcolor(gs12)) /// 
 (scatter empshare_prod_industry1 year, mcolor(gs4) msymbol(O) legend(label(3 "Below median"))) /// 
 (scatter empshare_prod_industry2 year, mcolor(gs12) msymbol(T) legend(label(4 "Above median "))), /// 
 xlabel(1880 (20) 2020) ytitle("Industrial employment share") legend(order(- "Industrial employment" "share 1882:" 4 3)) name(indempmed1882, replace) graphregion(color(white))
 
graph export "$reversing/results/figures/figA4.eps", replace
	
restore	




***************
*** Appendix Figure A-5: The effect of early industrialization on the number of universities and colleges, 1900-2019

** Panel A: All universities and colleges

** Choose covariates
local covariates log_land_access1 town_1700_perarea

local years "1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010"
matrix C90 = J(3,12,.)
matrix C95 = J(3,12,.)
local j = 1
foreach year in `years' {
		acreg unis_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100) 
		lincom empshare_ind_1882_std, level(90)
		matrix C90[1,`j']= r(estimate) \ r(lb) \ r(ub)
		lincom empshare_ind_1882_std, level(95)
		matrix C95[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}


local colnames
foreach year in `years' {
    local colnames `colnames' `year'
}
matrix colnames C90 = `colnames'
matrix list C90

matrix colnames C95 = `colnames'
matrix list C95

coefplot (matrix(C95), ci((C95[2] C95[3])) mcolor(black) ciopts(lpattern(solid) lcolor(black) msize(medlarge) recast(rcap))) (matrix(C90), ci((C90[2] C90[3])) mcolor(black) ciopts(lpattern(solid) lcolor(gray) lwidth(medium) mcolor(gray) msize(medlarge) recast(rcap))), xtitle("Year") yline(0, lcolor(black)) vertical graphregion(color(white)) at(_coef) xlabel(1900(10)2010, nogrid) ytitle("Number of large universities") ylabel(-0.60(0.20)0.80, nogrid) legend(off)

graph export "$reversing/results/figures/figA5a.eps", replace

** Panel B: Large universities and colleges

** Choose covariates
local covariates log_land_access1 town_1700_perarea

local years "1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010"
matrix C90 = J(3,12,.)
matrix C95 = J(3,12,.)
local j = 1
foreach year in `years' {
		acreg lunis_`year' (empshare_ind_1882_std = log_coal_access1) `covariates' if year == 2019, latitude(latitude) longitude(longitude) spatial bartlett distcutoff(100) 
		lincom empshare_ind_1882_std, level(90)
		matrix C90[1,`j']= r(estimate) \ r(lb) \ r(ub)
		lincom empshare_ind_1882_std, level(95)
		matrix C95[1,`j']= r(estimate) \ r(lb) \ r(ub)
		local ++j
}


local colnames
foreach year in `years' {
    local colnames `colnames' `year'
}
matrix colnames C90 = `colnames'
matrix list C90

matrix colnames C95 = `colnames'
matrix list C95

coefplot (matrix(C95), ci((C95[2] C95[3])) mcolor(black) ciopts(lpattern(solid) lcolor(black) msize(medlarge) recast(rcap))) (matrix(C90), ci((C90[2] C90[3])) mcolor(black) ciopts(lpattern(solid) lcolor(gray) lwidth(medium) mcolor(gray) msize(medlarge) recast(rcap))), xtitle("Year") yline(0, lcolor(black)) vertical graphregion(color(white)) at(_coef) xlabel(1900(10)2010, nogrid) ytitle("Number of universities") ylabel(-0.60(0.20)0.80, nogrid) legend(off)

graph export "$reversing/results/figures/figA5b.eps", replace

*** EOF

